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Using recent dimensionality reduction techniques in large systems of coupled phase oscillators exhibiting bistability, we analyze 
complex macroscopic behavior arising when the coupling between oscillators is allowed to evolve slowly as a function of either 
macroscopic or local system properties. For example, we observe macroscopic excitability and intermittent synchrony in a system 
of time-delayed Kuramoto oscillators with Hebbian and anti-Hebbian learning. We demonstrate the robustness of our findings by 
considering systems with increasing complexity, including time-delayed oscillators with adaptive network structure and community 
interaction, as well as a system with bimodally distributed frequencies. 



1. Introduction 

Large systems of coupled oscillators are ubiquitous in na- 
ture and serve as a basis to study collective behavior. Some 
examples include synchronized flashing of fireflies |1 1, cardiac 
pacemaker cells |2|, walker-induced oscillations of the Mil- 
lennium Bridge |3|, Josephson junction circuits lH, audiences 
clapping |5|, circadian rhythms in mammals f6\, cell function 
Q, neural processing |8|, and chemical oscillations ll9l [TOl . 
In certain situations these oscillators can be approximately de- 
scribed in terms of only their phase angle 9. Kuramoto showed 
im that the evolution of the phases in an ensemble of N weakly 
coupled oscillators approximately obeys 
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where 6'„ and aj„ are, respectively, the phase and intrinsic fre- 
quency of oscillator n, and Hnm is a 27r-periodic function that 
describes the coupling between oscillators n and to. When 
such oscillators represent limit cycles arising from a Hopf bi- 
furcation their coupling is generically sinusoidal, leading to 
the choice Hnm{S) — {knm / N) sm{9) . When coupling is 
uniform, i.e. knm = k, one obtains the classical Kuramoto 
model which has become a paradigm for the study of synchrony 
in coupled heterogeneous oscillators. Generalizations of the 
Kuramoto model have become an important area of recent re- 
search, including investigations of non-sinusoidal coupling 1 1 1 1, 
cluster synchrony |12|, the effects of network topology Ill3l 
[T4ll . non-local coupling [15], external forcing [161, coupled ex- 
citable oscillators IITTI . phase resetting ifTSl . time-dependent 
connectivity lfT9]| . noise ll20l . and communities of coupled os- 
cillators [2Tl|22l- Recently, the analysis of many such systems 
has been simplified by a dimensionality reduction proposed by 
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Ott and Antonsen ll23l l24l . making many cases analytically 
tractable for the first time. 

A major difficulty in the study of complex systems (e.g., 
neural processing and cell function) is overcoming the common 
disconnect between simple microscopic and complex macro- 
scopic dynamics referred to as emergence. In this paper, we 
study emergent macroscopic behavior that cannot be deduced 
from the individual oscillator dynamics alone, but requires a 
systems-level analysis. We study macroscopic dynamics arising 
when slow coupling adaptation is combined with large systems 
of oscillators. Natural examples of systems involving adaption 
of system parameters such as coupling strength include clap- 
ping audiences [Q, brain fluctuations ID, regulation of sleep 
and circadian rhythms f25l, and regulation of cardiac behav- 
ior [26] . One natural way to model adaptive dynamics in such 
systems is to allow for the connectivity to evolve as a func- 
tion of the degree of synchrony of the system. Recent studies 
on adaptive oscillator systems have largely modeled two types 
of synaptic plasticity: spike-timing-dependent plasticity ll27l 
and Hebbian learning [28]. We further classify such adapta- 
tion rules as either uniform adaptation (the evolution of global 
coupUng in all-to-all coupled systems according to global sys- 
tem properties) or network adaptation (the evolution of individ- 
ual links in possibly heterogeneous networks according to local 
properties), both of which are studied here. 

The inclusion of adaptive rules in Kuramoto-type systems 
can result in rich dynamics that has sometimes been proposed 
to model information processing in the brain !l27l |281 . Typi- 
cally, however, these adaptive rules are added to the standard 
Kuramoto model, which has relatively simple macroscopic dy- 
namics (e.g. no memory). In this paper, we explore the addi- 
tion of adaptive rules to oscillator systems exhibiting bistability, 
such as those studied in Refs. [29., 30, 3h 32} . When com- 
bined with adaptation, we find that bistability allows for com- 
plex macroscopic behavior such as excitable and intermittently 
synchronous states in addition to simple steady-state behavior. 
In this paper we consider the case where the timescale of cou- 
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pling adaptation is much larger than the timescale of oscillator 
dynamics, which will allow us to separate time scales, first solv- 
ing for the fast oscillator dynamics using the work of Ott and 
Antonsen |23| and then analyzing the slow adaptation dynam- 
ics. We find that even when the adaptation is chosen to be a 
simple function of the system state, a variety of macroscopic 
behaviors can be attained by varying parameters of that simple 
function. The dynamics described in this paper fall within the 
framework of dynamic bifurcation theory |33 1, which describes 
bifurcations that occur in fast dynamics in response to one or 
more slowly-changing parameters. In this paper the bifurca- 
tions correspond to transitions between macroscopic incoher- 
ent and synchronized states in response to one or more slowly 
changing coupling strengths. 

This paper is organized as follows. In Sec. 2 we study a 
system of time-delayed oscillators subject to uniform adapta- 
tion. In Sec. 3 we study three more complicated models that 
yield richer dynamics: (A) network adaptation on a system of 
time-delayed oscillators, (B) community interaction between 
two communities of time-delayed oscillators with community- 
wise uniform adaptation, and (C) uniform adaptation on a sys- 
tem of oscillators with intrinsic frequencies drawn from a bi- 
modal distribution. Finally, in Sec. 4 we conclude by discussing 
our results. 

2. Time-delayed oscillators with uniform adaptation 

In this section we study a system of N oscillators coupled 
through a time-delayed order parameter in which the coupling 
strength is allowed to slowly adapt in response to the values of 
this order parameter. This system allows for analytic results that 
will later serve as a guide to the analysis of more complicated 
systems. Letting aj„ denote the intrinsic frequency of oscillator 
n [which we assume to be randomly drawn from a distribution 
g{(jj)] and r ^ J2n=i denote the Kuramoto order pa- 
rameter, we consider the following model. 



On = ujn + klm{ze 
Tz = r — z. 



(2) 

(3) 



Since Eq. (3 i can be written as z{t) = r^^ r{t')e^^ dt' , 
z may be interpreted as a time-delayed version of r. This time- 
delayed order parameter is the one that affects individual oscil- 
lators. In the continuum limit, N ^ oo, this system has been 
shown to represent exactly the case where the coupling between 
pairs of oscillators in the Kuramoto model is time-delayed with 
time delays that have an exponential distribution with average 
r li29l . Note that t — > yields z = r, which recovers the 
Kuramoto model |9|. We extend this system by allowing the 
uniform coupling constant k to slowly adapt following 



Tk = G{k,z), 



(4) 



where T is the timescale of adaptation and G is a function that 
describes the adaptation of k in terms of its current value and 
the perceived (delayed) order parameter. We will assume that 
T is much larger than both r and the time scale of oscillator 



dynamics, given by the inverse of the spread of g{uj) ||24| so that 
we may utilize a separation of time scales to solve Eqs. Q and 
^ assuming constant k and then solve Eq. Q while assuming 
a steady state. Note that letting T — > cx) recovers the non- 
adaptive system with fixed k. 

2.1. Fast oscillators in the continuum limit 

We begin by solving Eqs. (j2| and ([3]l for fixed k in the con- 
tinuum limit. We let /(w, f?, t) denote the density of oscillators 
with frequency lo and phase 9 at time t. Conservation of oscilla- 
tors implies that /(w, 6*, t) must satisfy the continuity equation 



dtf + de{fe) = Q. 



(5) 



Following Ref. |!291, this partial differential equations (PDF) 
can be reduced to the single complex-valued ordinary differen- 
tial equation (ODE) 
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where we have assumed the frequency distribution is Lorentzian, 
i.e. g{Lo) = A7r^^/[A^ + (w — wq)'^]- This assumption is not 
necessary in general but simplifies the presentation of our re- 
sults. Along with Eq. ([3]), Eq. ^ closes the oscillator dynamics. 
We note that Eq. (j6]l was derived in Ref. L29il for time-delayed 
oscillators without coupling adaptation. 

Assuming a fixed k value, we now look for steady state so- 
lutions by defining r — Re^^, z = pe"^, and setting R — p — 
and ip — (/) — H.. Without loss of generality, by rescaling time 
we can set A 1. As shown in Fig.[T]for ui^ — 5, in addition 
to the incoherent solution R ^ p — 0, a pair of synchronized 
solutions appears at fci = 2wo, given by ll34ll 
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with a corresponding angular velocity given by 
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Subscripts sju denote whether the solution is stable or unsta- 
ble, respectively. AtA;2 = (wq + 4) /2 the unstable synchronized 
branch merges with the incoherent solution, which becomes un- 
stable for k > k2- Note that for k between ki and /c2 we find 
bistability since there are both coherent and incoherent solu- 
tions that are stable to perturbation (the linear stability of these 
solutions has been discussed in Refs. ||29l [34l ). Furthermore, 
along the synchronized branches cf) lags behind t/; by an angle 
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We now remark on two aspects of this system not previ- 
ously discussed in Ref. 1291 1341 . First, we note that the average 
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Figure 1 : (Color online) Solutions Rg (solid blue curve) and _R„ (dashed red 
curve) [given by Eq. ^] for the time-delayed system with ojq = 5 and A = 1. 
Inset: Sis (solid blue curve) and Qu (dashed red curve) corresponding to the 
angular velocities of the synchronized states. Note that Qs and Qu are much 
smaller than the average intrinsic frequency uiq = 5 (dotted line). 



angular velocity il of oscillators in a synchronized state is con- 
siderably less than the average intrinsic frequency ujq (see inset 
in Fig. [T]l. This reflects the fact that each individual oscillator 
is coupled not to the instantaneous mean-field, but to the time- 
delayed version, slowing down the entire synchronized popula- 
tion. 

Second, whereas the distribution of locked oscillators in the 
standard Kuramoto model is symmetric about the mean oscil- 
lator frequency ujq, this symmetry is broken by the time de- 
lays and as a result the distribution of locked oscillators for 
oscillators with time-delayed coupling is biased toward oscil- 
lators with angular frequencies near 51. Because fl is much 
smaller than ojq, this distribution of locked frequencies is typi- 
cally spread asymmetrically around the mean frequency uq. We 
compute the critical frequencies ujc,± separating phase-locked 
and drifting oscillators by entering a rotating frame in which 
synchronized oscillators appear stationary by defining 9„ ~ 
On — 4>- Here 0„ evolves according to 0„ = cj„ — — 
fcpsin(0„), so that 0„ reaches an equilibrium and becomes 
phase-locked if |aj„ — fi| < kp, and otherwise drifts indefi- 
nitely. Thus, the critical frequencies that separate the drifting 
and locked populations are ujc.± = 51 ± kp. 

2.2. Slow coupling adaptation 

Having solved the oscillator dynamics that evolve on the 
fast time scale, we now study adaptation given by Eq. (|4]) that 
evolves on a slow time scale. For simplicity, we assume that k 
relaxes to a linear function of p, 



G{k,p) ^a + /3p-k, 



(11) 



and will study the resulting behavior as a function of the pa- 
rameters a and /3. While this form for G is not essential, it sim- 
plifies our exploration of complex macroscopic behavior under 
adaptive uniform coupling while yielding rich dynamics. Using 
Eqs. Q and (|8]l, the behavior of the order parameter magnitude 
p and coupling strength k is described on the slow time scale 
by the ODE 



when the system is synchronized, and 

Tk ^ a — k. 



(13) 



when the system is incoherent (i.e., p — 0). 

As shown in Fig. |2|a), when the system is in the incoher- 
ent state and k surpasses fc2 a dynamic bifurcation occurs in 
a rapid transition from incoherence to synchronization. Simi- 
larly, when the system is in the synchronized state and k de- 
creases below ki another dynamic bifurcation occurs in a rapid 
transition from synchronization to incoherence. These rapid 
transitions from one branch to the other represent discontinuous 
phase transitions, which have also been refered to as explosive 
synchronization |[35l . Again, T is assumed to be large enough 
that state transitions occur with fixed k. Furthermore, we as- 
sume that, upon a perturbation of the oscillator phases which 
changes the value of p and R, the system returns to the inertial 
manifold in which the Ott-Antonsen ansatz is valid on a time 
scale which is much faster than T, so that we can assume k is 
constant during this process. It follows that the macroscopic 
behavior depends on the location of the stable fixed points of 
Eqs. ^V2\ and (13 1 as depicted for various situations in Figs.[2|a)- 

We now classify the nature of macroscopic behavior by study- 



ing the stable fixed points of Eq. ( 12 1 (the synchronized fixed 



point, fc* ) and Eq. ( 13 i (the incoherent fixed point, fc*„^). 
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Note that incoherent fixea points for k > k2 are not relevant 
since the incoherent branch is unstable in that region, and there- 
fore we will ignore these fixed points in what follows. Classify- 
ing the macroscopic dynamics for a particular choice of (a, /3) 
reduces to the analysis of the stable fixed points of Eqs. ( [T2) l 
and (13 I subject to the incoherent solution {p = 0) and stable 
synchronized solution [Eqs. ([8])]. The different possible macro- 
scopic behaviors are the following: 

• No synchronized fixed point, no incolierent fixed point. 

As shown in Fig. |2|a), if no stable fixed point exists on 
either branch the system will repeat the following macro- 
scopic oscillation: k will increase along the incoherent 
branch, then after the system synchronizes at k2, it will 
decrease along the synchronized branch, until desynchro- 
nization occurs at ki . We define this behavior as the in- 
termittent state and investigate its properties in Sec. IIC. 

• Synclironized fixed point, no incolierent fixed point. 

If a stable fixed point occurs only on the synchronized 
branch, we define two subclasses for this state: if A:*^^^^ > 
k2 we define the macroscopic state as the synchronized 
state [see Fig. |2j^b)], whereas if < fc2 we refer 

to the state as the excitable synchronized (ES) state [see 
Fig. |2jc)]. This distinction is made to account for the 
possibility that perturbations to the value of p in the ES 
state (e.g., due to noise or finite-size fluctuations [361 ) 
may desynchronize the system by decreasing p below pu, 
which would result in p — J- followed by k increasing 
until k = k2, after which the system will synchronize and 
return to the fixed point. Thus, the synchronized state can 
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Figure 2: (Color online) Various macroscopic behaviors may occur for unifoiTn adaptation following Eq. \\\\ depending on a, /9 and the location of fixed points 
(circles). Examples include (a) intermittent, (b) synchronized, (c) ES, and (d) bistable states. 



be interpreted as the resting state of an excitable system 
and a temporary desynchronization as an excitation. 

• Incoherent fixed point, no synclironized fixed point. 

Analogous to the previous case, we define two sub-classes 
for this state: if k*^^^ < ki we refer to this state as the in- 
coherent state, whereas if > ki we refer to this state 
as the excitable incoherent (EI) state. Again this distinc- 
tion is made to account for the possibility of perturba- 
tions to the value of p in the EI state that can produce a 
temporary synchronization. In this case the system can 
be synchronized if p is increased above p„, resulting in 
p —> ps followed by k decreasing until k = ki, desyn- 
chronization, and finally a retum to the fixed point. In this 
scenario, the incoherent fixed point can be interpreted as 
the resting state of an excitable system and a temporary 
synchronization as an excitation. 

• Synclironized fixed point, incolierent fixed point. If 

stable fixed points occur on both branches, we refer to 
this state as the bistable state, an example of which is 
shown in Fig.[2|d). 

To find the location of the bifurcations between these states, 
we calculate the critical a, (3 that correspond to the formation or 
destruction of fixed points on either branch. For fixed points on 
the incoherent branch, this occurs at a — k2- For fixed points 
on the synchronized branch, we require that the curves k — a 
and (3ps{k) are tangent if /3 > 0, which occurs when 



,{k) = k — a, 



, dpa{k) 
dk 



1, 



and coincide if /3 < 0, which happens when 

a + I3ps{ki) = ki. 



(14) 



(15) 



Finally, the boundary between EI and incoherent states is given 
by a = fci (the incoherent fixed point entering the bistable re- 
gion), while the boundary between ES and synchronized states 
is given by the curve a + /3ps (^2) = fc2 (the synchronized fixed 
point entering the bistable region). In Fig. |3] we show the bi- 
furcation diagram for loq — b and A = 1 by plotting curves 
describing the formation/destruction of incoherent fixed points 
in solid blue, synchronized fixed points in dashed red, and the 



40 



20 



-20 



-40 



Bistable 








Synchronized 




■^^\ 


\ ES 


Incoherent 








EI 








Intermittent ^^^^^ 



10 15 20 25 30 



a 

Figure 3: (Color online) Biforcation diagram summarizing boundaries between 
intermittent, synchronized, ES, incoherent, EI, and bistable states for ujq = 5 
and A = 1. 



borders between EI/ES and Incoherent/Synchronized states in 
dotted black. We label regions with the states described above. 
We note that excitable and intermittent states are possible only 
when /? < 0, which we refer to as anti-Hebbian adaptation 
(accordingly, we refer to > as Hebbian adaptation). This 
terminology is based on the observation that for /3 > (/3 < 0) 
in Eq. (Ill, coupling is promoted (inhibited) by the synchrony 



of oscillators. 

2.3. Intermittent case 

Due observations of intermittent synchrony in various os- 
cillator systems (e.g., in neural activity 1 8 , 37 1 and clapping au- 
diences ||5l) we now study in detail the intermittent case illus- 
trated in Fig. [2]^a) and characterized by intermittent periods of 
macroscopic synchronization. Of interest is the period of oscil- 
lation, which can be found by integrating the time spent follow- 
ing the incoherent and synchronized branches of the bistable re- 
gion. The time spent in the incoherent state. Tine, corresponds 
to the time it takes for k to increase from ki to ^2 with p = 
and is given by 
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Figure 4: (Color online) (a) Theoretical Ti^op (solid green), Ti„t, (dashed red), and Tsync (dot-dashed blue) agree well with Ti^op averaged over 16 simulations 
(black asterisks, where error bars indicate standard deviation). Other parameters are (3 = —20, a;o = 5, A = 1, T = 2000, and A'^ = 5000. (b) For a = 18 
theoretical values for R{t) (solid blue) agree well with direct simulation \r{t)\ (red crosses). 



Similarly, the time spent in the synchronized state, Tsync, cor- 
responds to the time it takes for k to decrease from k2 to ki 
along the synchronized branch and is given by 



T —T 

sync ^ 



a 



dk 



(17) 



Since we assume that the timescale of adaptation is much larger 
than the timescale of oscillator dynamics, we neglect the time it 
takes for oscillators to synchronize and desynchronize at fc2 and 
fci, respectively. This gives the period of oscillation Tioop ~ 
Tsync + T^nc- Fig. a) shows Tioop (soHd green curve). Tine 
(dashed red curve), and Tsync (dot-dashed blue curve) as a func- 
tion of a for Wo = 5, A = 1, and /3 = —20. For these param- 
eters fci = 10 and k2 — 14.5. In addition, we compute the pe- 
riod of oscillation from simulating N = 5000 oscillators with 
T — 2000, plotting the mean of Tioop over 16 simulations at 
each a (black asterisks). Error bars indicate the standard devia- 
tion. While Tine and Tioop diverge d& a ^ k^ , Tsync and Tioop 
remain finite as a — > [fci — j3ps(ki)\'^, since the square root 
singularity of Ps(fc) at fc = fci prevents the integral in Eq. ( 17 1 
from diverging. 

As shown in Fig. Qb), the macroscopic behavior of the 
system oscillating between incoherent and synchronized states 
may be described by considering the low dimensional system 
given by Eqs. (j7|i, (|8]l, (12i, and (13i. This theoretical solu- 
tion R{t) (solid blue curve) agrees well with the order param- 
eter's magnitude \r{t)\ (red crosses) from direct simulation of 
the high-dimensional system given by Eqs. (|2]i, ([3]), (|4]), and 



(111. The simulation in Fig. |4|b) was done with N = 5000 os- 
cillators with T = 2000, a = 18, and /3 = -20. Remarkably, 
the behavior of the high-dimensional system is captured well 
by this piecewise defined one-dimensional ODE. The period 
taken from simulations is slightly longer than our theoretical 
solution, which is most likely due to two effects. First, our the- 
oretical solution neglects the synchronization and desynchro- 
nization times at the dynamic bifurcations occurring at fc = fci 
and fc2. Second, along the incoherent branch the value of the 
order parameter in simulations typically takes values of size 
O (AT- 1/2) [36J rather than zero, which slightly slows down the 
adaptation. 
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Figure 5: (Color online) Spiking events are shown for the order parameter r 
when a = 14.45 and /3 = —20 are chosen so that our system is in the EI 
state. Note that the timescale of spiking is dominated by the spontaneous syn- 
chronization process |31| . (inset) For a system of size N we may predict the 
expected time between spikes as T^pike °^ exp(f A^) for some constant C 
The solid line indicates a least-squares fit Tspike 

IT = 2.16 exp(0. 00217V). 



2.4. Excitable incoherent case 

We conclude our analysis of this model by studying the EI 
state. As previously mentioned, if the system is in the incoher- 
ent state with fc = ^Xnc^ ^ perturbation to the order parameter 
can cause r to become larger than the unstable solution r„, re- 
sulting in a dynamic bifurcation. While fc remains fixed during 
this rapid transition, after synchrony fc will evolve until the sys- 
tem returns to the equilibrium of r ~ and fc = fc*„c- In 
particular, for finite systems this perturbation could occur due 
to finite size effects, resulting in a spontaneous synchronization 
event |31|. This can be viewed as a random spiking event for 
the macroscopic dynamics, which corresponds to the oscillators 
synchronizing very briefly relative to the typical time between 
spikes. 

Spiking events are shown in Fig. [S] where we plot |r(i)| 
versus time for r = 1, T = 1000, A = 1, = 5, and 
TV — 1050. Note that the system spends the majority of time 
in the incoherent state, and the slow timescale of spontaneous 
synchronization dominates other time scales. Defining the av- 
erage time between synchronization events as the inter-spike 
time, Tspike, we briefly discuss the dependence of Tspike on 
system size iV. In Ref. ||3T1 it was shown that the sponta- 
neous synchronization event can be modeled as a Kramer es- 
cape process where the expected escape time is proportional 
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Figure 6: (Color online) Example {k, \r\) trajectories of the system given by 
Eqs. ^19|-j2TJ. The solid blue and dashed red trajectories were obtained using 
{a, p) = (24, 0) and (6, 0), respectively with an initial coupling strengths 
fc = 6 and 22, respectively. Other parameters are = 5, A = 1, t = 1, 
T = 1000, 7 = 3, dmin = 100, and = 1000. 



to exp(CA^) for some constant (. Therefore, because the es- 
cape process dominates the timescale of dynamics, we expect 
that the inter-spike time scales as Tgpike exp(CA^). This 
is confirmed in the inset of Fig. |5j where Tgpike is shown to 
vary exponentially with N. The solid line is a least-squares fit 

3. Other models 

In the previous section we analyzed in detail the model 
given by Eqs. Q-Q, which describes a system of oscillators 
with heterogeneous time-delays subject to uniform coupling adap- 
tation. The purpose of this model was to illustrate generic be- 
havior occuring in adaptive networks with bistable regimes. In 
this section we study numerically and analytically several other 
models which have been selected to show that the type of be- 
haviors observed in the previous section occur more generally. 
In particular, in Sec. 3.1 we investigate network adaptation, 
which is often used in Kuramoto-type models of information 
processing and memory in neural networks Il27ll28]| . In Sec. 3.2 
we explore complex macroscopic behavior that can arise for 
adaptation in networks containing community structure. Fi- 
nally, in Sec. 3.3 we show that our findings apply to other os- 
cillator systems exhibiting multistability (e.g., due to frequency 
adaptation [311 or inertia |32|) by studying adaptation in oscil- 
lator systems with bistability due to a bimodal distribution of 
intrinsic frequencies |30|. 

3.1. Network adaptation 

First, we will consider a system similar to Eqs. Q-Q in 
which the interactions between oscillators are not mediated by 
a global mean field, but occur instead through an underlying 
coupling network. We assume the undirected network structure 
is represented by an adjacency matrix A, where 



A — 



1 if a link exists from oscillator m to oscillator n, 
if no link exists. 
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Figure 7: (Color online) Bifurcation diagram summarizing oscillatory (black 
squares), synchronized (red circles), ES (cyan asterisks), incoherent (blue tri- 
angles), EI (green plusses), and bistable (yellow crosses) states for network 
adaptation of time-delayed oscillators with A = 1, ujq = 5, t = I, and 
T = 1000. 



Introducing a coupUng weight fc„„j to each link and using the 
locally-defined order parameters r„, where 



N 



we consider the system given by 

9n = + A^^Im(z„e"'''"), 



T Zj^^ Tn^l Z^y-, 



Tkri 



a H- /3Re{rnZ^) - 



(18) 



(19) 
(20) 

(21) 



where a;„ is again randomly drawn from a Lorentzian with 
mean luq and spread A, and Ad is the dominant eigenvalue 
of A. We normalize the coupling term in Eq. (19 1 by Xjj so 
that the /(;„„, values producing bistability are on the same or- 
der as k values that yield bistability in the uniform adaptation 
model Il3l . To measure the global degree of synchrony and 
coupling strength we introduce the average order parameter 



V A 
and average coupling strength 



e [0,1], 



nm^nm 



(22) 



(23) 



Since we are assuming that oscillator n is affected by a de- 
layed order parameter, the adaptation of the coupling fc„,„ be- 
tween oscillators m and n, Eq. (21 1, is assumed to depend on 
the local instantaneous order parameter of oscillator n, rn, and 
the delayed order parameter at oscillator to, z„j. As before, we 
interpret positive values of /3 as Hebbian adaptation, and nega- 
tive values as anti-Hebbian adaptation. 

Using the Chung-Lu model I.38J . we construct an undirected 
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Figure 8: (Color online) Community interaction model with parameters A = 1, ujq = 5, a = 18, /9 = —30, T = 200, and e = 0.105 (a) and 0.085 (b). Top 
panels: evolution of |ri | (solid blue line) and |r2 1 (dashed red line), bottom panels: evolution of fcn (solid blue line), ^22 (dashed red line), fci2 (dot-dashed green 
line), and k2i (dotted black line). 



network with a power-law degree distribution, P{d) oc d^"', 
with exponent 7 = 3, minimum degree d„iin — 100, and 
N ~ 1000 oscillators, where the degree d is defined as d„ = 
X]m=i ^nm- The parameters for the oscillator dynamics are 
Wo = 5, A = 1, T = 1, and the adaptive timescale is T = 1000. 
The dominant eigenvalue for the network constructed for the 
simulations shown here is Ad — 232.325. In Fig.|6]we show 
representative (fc, |f|) trajectories. First using (a,/?) = (24,0), 
we allow the average coupling strength to increase from an ini- 
tial value of fc = 6 (solid blue trajectory). Next using (a, (3) — 
(6, 0), we allow k to decrease from an initial value of /c = 22 
(dashed red). We find that in analogy with the uniform adapta- 
tion case, a stable synchronized solution \R\ > is created at 
k = ki ~ 12.6, and the incoherent solution |_R| = becomes 
unstable at /c = ^2 ~ 13.6. Thus, dynamic bifurcations occur 
approximately when an incoherent state's average coupling in- 
creases through k2 or a synchronized state's average coupling 
decreases through fci. 

Next we numerically explore the {a, (3) parameter space, 
classifying the observed behaviors as bistable, intermittent, syn- 
chronized, ES, incoherent, and EI, following the criteria in Sec. 
2. In Fig.|7]we plot the results. Starting from the top left and 
proceeding clock-wise, we plot bistable (yellow crosses), syn- 
chronized (red circles), ES (cyan asterisks), oscillatory (black 
squares), EI (green plusses), and incoherent (blue triangles). 
These states were found by tracking the trajectories of |f | and k 
for two simulations at each pair [a, (5), one trajectory starting 
from an incoherent state with k < ki, and the other starting 
from a synchronized state with k > k2- The results are smooth 
enough so that boundaries between regions are clear As ex- 
pected, while the exact boundaries in Fig. |7] differ from those 
plotted in Fig. |3] the topologies of the two phase spaces agree 
qualitatively. 

3.2. Community interaction 

Next, we generalize the system studied in Sec. 2 to a two 
community model where coupling is strong within communi- 
ties and weak between communities. For simplicity, we assume 
adaptation within and between each community is uniform. The 



model we consider is: 



Tk,^, ^G''^'{k,f,z), 



a' = l 



(24) 

(25) 
(26) 



where = 1, 2 denotes the community, denotes the phase 
of an oscillator in commuiuty a, Vc — X]m=i ^^^^ '■^^ 
Kuramoto order parameter over oscillators in community cr, 
k = [kii,ki2,k2i,k22]'^,r= [ri,r2]'^, z= [2:1,2:2]^, and the 
natural frequencies are drawn from the distribution .9^(0;). 

Separating the fast oscillator dynamics from the slow adap- 
tation dynamics as before, a dimensionality reduction for the 
iV^ ^ 00 limit as in Refs. ||2l]|23] yields 



1 ^ 

a' = l 



ka^:{z„, - zl,rl), (27) 



where we have assumed that the distribution {lu) is Lorentzian 



with spread and mean Wq . Eqs. (27 1 and (25 1 give the low- 
dimensional evolution of oscillator dynamics. Furthermore, we 
consider the adaptation dynamics given by 



(28) 



Depending on the choices of A^, ujq, aa-a', and f3cra', the 
resulting dynamics can vary greatly. For simplicity we choose 
Ao- = A = 1 and Uq ^ utQ ^ 5, ^ t ^ 1, a^a' ~ ct and 
/3cra' = P for cr = cr', and aaa' — ^ct and j3aa' = e/3 for a ^ a' 
where < e < 1. We induce oscillatory behavior by choosing 
= 18, /3 = -30, and T = 200, and investigate the effect of 
varying e. Particularly, we are interested in macroscopic syn- 
chrony of the two communities. 

We simulate the system with — 2000 oscillators in 
both communities with initial coupling strengths of fcn = 9, 
^22 = 14, and kyi = k2\ — for values e — 0.105 and 
0.085. In Fig.lsla) we plot |ri(i)| (solid blue curve) and \r2(t)\ 
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Figure 9: (Color online) (a) Bifurcation diagram for the Kuramoto model with bimodal frequency distribution. Transcritical, Hopf, homoclinic, and saddle 
node/SNIPER bifurcations are plotted in dashed back, clue, green, and red, respectively, paths taken in Fig. |10| are ploted in solid black, (b) Zoomed-in view 
of bistable regions. We label regions where incoherent, synchronized, and standing-wave solutions are stable In, S, and SW, respectively. 



(dashed red curve) in the top panel and fcn (t) (solid blue curve), 
^22 (t) (dashed red curve), ki2 {t) (dot-dashed green curve), and 
k2i{t) (dotted black curve) in the bottom panel for e — 0.105. 
In Fig. |8jb) we plot the same quantities for e = 0.085. Al- 
though the two communities start in out-of-phase macroscopic 
states, for e = 0.105 the macroscopic dynamics of the commu- 
nities synchronize near t — 3T. However, for e = 0.085 they 
remain out of phase past t ~ 200T. This significant difference 
in behavior for such a small change in e suggests a sensitive 
dependence on the system parameters in addition to initial con- 
ditions. 

3.3. Bimodal frequency distribution 

Finally, we study uniform adaptation of a system of oscilla- 
tors without time-delay, but having bistability due to a bimodal 
distribution of intrinsic frequencies. The model we study is the 
following: 



An = w„ + fclm(re *^"), 
Tk = G{k,z). 



(29) 
(30) 



where r — En=i normal Kuramoto order pa- 

rameter and now we assume w„ are drawn from the double 
Lorentzian 



A 

2^ 



(w-wo)2 + A2 (w + CJo)2 + A2 



(31) 



which is bimodal for A < \/3t^o- We note that in Ref. 
a similar oscillator system with bimodally-distributed frequen- 
cies is studied, but with an explicitly time-dependent sinusoidal 
couphng strength rather than system-dependent coupling adap- 
tation. 

This model is particularly interesting because in addition to 
the simple coherent and incoherent fixed points, stable solutions 
can also take the form of standing waves in which two synchro- 
nized groups [one corresponding to each peak of g(w)] oscillate 
with opposite angular velocity ll30l . These solutions are found 
for intermediate coupling strengths such that groups of oscilla- 
tors with frequencies near ujq and —utQ synchronize, but these 



two groups do not synchronize with one another These two 
groups act as giant oscillators that continue to pass one another, 
maximizing | r \ when the two groups have equal phase and min- 
imizing |r| when they have opposite phase. For a detailed anal- 
ysis of this oscillator dynamics refer to Ref. |30|. 

In Fig. |9|a) we summarize the bifurcation diagram. Hor- 
izontal and vertical axes are Aojo/k and 4A/fc, respectively, 
and transcritical, Hopf, homoclinic, and saddle-node/SNIPER 
bifurcations are plotted in dashed black, blue, green, and red 
curves, and labelled TC, HB, HC, and SN/SNIPER, respec- 
tively. In Fig. [9|b) we show a zoomed-in view of the bistable 
regime and indicate regions where the incoherent, synchronized, 
and standing-wave solutions are stable. Regions are labelled S, 
In, and/or SW if the synchronized, incoherent, and/or stand- 
ing wave solutions are stable in that region, respectively. For 
small K the incoherent solution is the only stable solution. This 
solution loses stability either in a transcritical bifurcation or a 
Hopf bifurcation, giving rise to synchronized or standing-wave 
solutions. Synchronized solutions are also born at the saddle- 
node/SNIPER bifurcations, and the standing-wave solution dis- 
appears at the homoclinic bifurcation. There are two distinct 
regions of bistability in the approximately triangular area in the 
middle of the plot. For 4A/A: > 1 [labeled S/In in Fig.|9];b)] the 
synchronized and incoherent solutions are both stable, whereas 
for 4A/fc < 1 [labeled S/SW in Fig.|9]:b)] the synchronized and 
standing- wave solutions are stable. 

Letting Eq. ( 30 1 take the linear form given by Eq. (Ill with 
a = 5, /3 = —5, and r = 1000, we simulate a system with N = 
2000 oscillators for wo = 1 and (a) A = 0.82, (b) 0.89, and (c) 
1.02. The respective trajectories in phase space (see solid black 
lines in Fig.|9|l yield the following behaviors: (a) the system os- 
cillates between synchronized and standing wave states; (b) the 
system repeats a synchronized^- incoherent —>standing wave — >^ 
synchronized cycle; and (c) the system oscillates between syn- 
chronized and incoherent states. We plot the behavior of each in 
(fc, |r|) space in Figs. 10 a), (b), and (c). Note in Fig. 10 b) that 



the macroscopic dynamics transition from incoherent to stand- 
ing wave at fc « 3.65 (see arrow) as predicted by the Hopf 
bifurcation in Fig. |9] In this case we see three dynamical bi- 
furcations in the transitions from incoherent— >standing-wave, 
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Figure 10: (Color online) \r{t)\ vs k{t) trajectories for the Kuramoto model with bimodal frequency distribution with uniform adaptation following Eq. jTTJ for 
JV = 2000 oscillators, ojq = 1, a = 5, /3 = —5, r = 1000, and A = 0.82 (a), 0.89 (b), and 1.02 (c). A transition from incoherence to a standing wave solution 
in (b) is indicated by an arrow. 



standing-wave—> synchronized, and synchronized— >^incoherent 
states. 

4. Discussion 

We have investigated analytically and numerically the effect 
of slow coupling adaptation on models of coupled phase oscil- 
lators exhibiting bistability and characterized complex macro- 
scopic behavior that extends to other bistable phase oscillator 
systems where bistability arises (e.g., due to frequency adap- 
tation BTl or inertial terms lf32l ). In addition to states with 
simple macroscopic fixed points, we have observed for uniform 
coupling adaptation on bistable systems macroscopic excitable 
and intermittently synchronous states. We leave open the ex- 
ploration of further dynamics that may occur for systems ex- 
hibiting multi-stability. 

Besides considering only uniform coupling adaptation (i.e., 
allowing the global coupling strength of an all-to-all system to 
evolve depending on macroscopic system properties), we have 
also addressed network adaptation (i.e., allowing the links be- 
tween individual oscillators to evolve according to their local 
properties). Network adaptation allows for heterogeneities in 
evolving networks to be accentuated and is often more realis- 
tic (e.g., Hebbian learning in neural systems [28 1). However, 
we have found that even when the underlying network struc- 
ture is heterogeneous, which in turn promote heterogeneities in 
the coupling between oscillators, qualitatively similar macro- 
scopic behavior emerges, i.e. fixed points, excitable, and inter- 
mittently synchronous states. Although our results for this case 
are purely numerical, we note that our results from the uniform 
adaptation model describe more heterogeneous networks with 
network adaptation very well. The development of more ad- 
vanced methods for dimension reduction for heterogeneous os- 
cillator networks is an open area of research, although progress 
continues |40|. 

We also have considered uniform adaptation for systems 
with either community interaction or bimodal frequency dis- 
tributions. In the community interaction model we have found 
complicated behavior even for simple parameter assumptions. 
We hypothesize that changing the manner in which communi- 
ties interact and/or increasing the number of communities could 
lead to richer, more complicated dynamics, including chaotic 
macroscopic states. In the bimodal frequency distribution model. 



we have demonstrated new dynamic bifurcations corresponding 
to the transitions between standing-wave solutions and the typ- 
ical incoherent and synchronized states. 

This work also provides a strategy for reconciling the com- 
mon disconnect between microscopic behavior (i.e. individual 
oscillator dynamics) and macroscopic phenomena. In the sys- 
tems studied in this paper we have shown that entire popula- 
tions of oscillators can combine into a single functional unit. 
For example, a wide range of parameters yields intermittent 
synchronous dynamics, which we liken to clock-like behavior. 
Similarly, we liken the dynamics of excitable and bistable states 
to neuron-like firing and switch-like behavior, respectively. One 
interesting direction of future research motivated by the work 
presented in this paper is the study of even more complex sys- 
tems that are composed of many functional units in a hierarchi- 
cal organization. In particular, one could study systems built out 
of different kinds of functional units, for instance to understand 
the resulting dynamics when networks of clocks, neurons, and 
switches interact. Because of their analytic tractability and sim- 
plicity, we believe that the results presented in this paper could 
prove a useful tool for understanding the generic behavior of 
these complex systems. 
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